Smoothing models

Get data - smooth disturbed and control plots together - AGB gain

library(mgcv); library(data.table); library(gratia); library(marginaleffects); library(ggplot2); library(dplyr); library(tidyr); library(zoo); library(slider); library(here)

#####test Paracou####

obs <- read.csv("G:/My Drive/cesab bioforest/Data/aggregated_data_v9.csv")

setDT(obs)

obs=obs[obs$Site=="Paracou"]


setnames(obs, "Plot", "subplot")

obs$plot=gsub("_.", "",obs$subplot )

obs <- obs[!Year%in%1984:1990]

obs <- dcast(
  obs,
  subplot + plot + Year ~ variable,
  value.var = "value"
)


##### upload data 
load("../clim_by_site.rda")

clim_all=rbindlist(lapply(clim_by_site[["Paracou"]], function(x) x$climate[[1]]),
                   idcol = "plot")

#"Paracou", "Mbaiki", "Corinto", "SUAS", "Lesong", "Ulu Muda", "Sungai Lalang", "Tene", "Jari"

baseline <- clim_all[year >= 1981 & year <= 2010]


clim_stats <- baseline[, .(
  mean_tmax = mean(tmax),
  sd_tmax   = sd(tmax),
  mean_vpd   = mean(vpd),
  sd_vpd     = sd(vpd),
  mean_srad  = mean(srad),
  sd_srad    = sd(srad),
  mean_def  = mean(def),
  sd_def    = sd(def)
), by = month]

clim_all <- merge(clim_all, clim_stats, by = "month", all.x = TRUE)

clim_all[, `:=`(
  z_tmax = (tmax - mean_tmax) / sd_tmax,
  z_vpd  = (vpd  - mean_vpd)  / sd_vpd,
  z_srad = (srad - mean_srad) / sd_srad,
  z_def  = (def  - mean_def)  / sd_def
)]

#climate anomalies follow inventories. For Paracou, mean per year before 1996, two years after

clim_all[, year_bin :=
           ifelse(year <= 1995,
                  year,                     # keep original year
                  1997 + ((year - 1996) %/% 2) * 2   # 2-year bins
           )]


census_anom <- clim_all[,.(
  mean_z_tmax = mean(z_tmax, na.rm = TRUE),
  mean_z_vpd = mean(z_vpd, na.rm = TRUE),
  mean_z_srad = mean(z_srad, na.rm = TRUE),
  mean_z_def = mean(z_def, na.rm = TRUE)
),
by = .(plot, year_bin)]

#AGB values before 1991 are not consistent
census_anom=census_anom[year_bin>=1991] 


merged <- merge(
  obs,
  census_anom,
  by.x = c("Year","plot"), by.y=c("year_bin","plot")
)


control=c(1,6,11,13,14,15)

merged$Treatment <- ifelse(merged$plot %in% control, "control", "disturbed")

merged$Treatment <- factor(merged$Treatment)
merged$plot <- factor(merged$plot)
merged$subplot <- factor(merged$subplot)

merged=merged[, agb_gain:=agb_growth+agb_recr]

GAM k=5, k=8 and bs = “fs”

NAs for agb_growth - plots 13,14 and 15 miss data (?)

merged <- merged %>%
  group_by(subplot) %>%
  mutate(
    fitted_GAM5 = predict(
      gam(agb_gain ~ s(Year, k = 5), data = cur_data()),
      newdata = cur_data()
    )
  ) %>%
  ungroup()


p1=ggplot(merged, aes(x = Year, y = fitted_GAM5, color = subplot)) +
  geom_line(size = 1) +  # predicted trajectory
  geom_point(data = merged, 
             aes(x = Year, y = agb_gain, color = subplot),
             inherit.aes = FALSE, size = 2) +  # observed points
  theme_bw() +
  labs(x = "Year", 
       y = "agb_gain", title="GAM-k=5") + theme(legend.position = "none")


plot(merged$agb_gain, merged$fitted_GAM5)

p1.1 <- ggplot(merged, aes(Year, agb_gain - fitted_GAM5)) +
  geom_point() +
  ggtitle("GAM residuals")

p1 + p1.1

ggplot(merged, aes(mean_z_tmax, agb_gain - fitted_GAM5)) +
  geom_point() +
  geom_smooth(method = "lm", color = "red", se = FALSE) +
  facet_wrap(~ Treatment, scales = "free_y")+
  ggtitle("GAM5 residuals vs Tmax")

ggplot(merged, aes(x = Year, y = fitted_GAM5, color = subplot)) +
  geom_line(size = 1) +  # predicted trajectory
  geom_point(data = merged, 
             aes(x = Year, y = agb_gain, color = subplot),
             inherit.aes = FALSE, size = 2) +  # observed points
  facet_wrap(~ subplot, scales = "free_y")+
  theme_bw() +
  labs(x = "Year", 
       y = "agb_gain", title="GAM-k=5") + theme(legend.position = "none")

ggplot(merged, aes(mean_z_tmax, agb_gain - fitted_GAM5)) +
  geom_point() +
  geom_smooth(method = "lm", color = "red", se = FALSE) +
  facet_wrap(~ Treatment, scales = "free_y")+
  ggtitle("GAM5 residuals vs Tmax")

GAM k=8

merged <- merged %>%
  group_by(subplot) %>%
  mutate(
    fitted_GAM8 = predict(
      gam(agb_gain ~ s(Year, k = 8), data = cur_data()),
      newdata = cur_data()
    )
  )%>%
  ungroup()

p11=ggplot(merged, aes(x = Year, y = fitted_GAM8, color = subplot)) +
  geom_line(size = 1) +  # predicted trajectory
  geom_point(data = merged, 
             aes(x = Year, y = agb_growth, color = subplot),
             inherit.aes = FALSE, size = 2) +  # observed points
  facet_wrap(~ subplot, scales = "free_y")+
  
  theme_bw() +
  labs(x = "Year", 
       y = "agb_gain", title="GAM-k=8") + theme(legend.position = "none")

p11

plot(merged$agb_gain, merged$fitted_GAM8)

p11.1 <- ggplot(merged, aes(Year, agb_gain - fitted_GAM8)) +
  geom_point() +
  ggtitle("GAM residuals")

 
p11.1

ggplot(merged, aes(mean_z_tmax, agb_gain - fitted_GAM8)) +
  geom_point() +
  geom_smooth(method = "lm", color = "red", se = FALSE) +
  facet_wrap(~ Treatment, scales = "free_y")+
  ggtitle("GAM8 residuals vs Tmax")

GAM5 bs = “fs”

#test GAM that is done on the full models (bs = "fs")

merged <- merged %>%
  mutate(
    fitted_GAM5_fs = predict(
      gam(agb_gain ~ s(Year, subplot, bs = "fs", k = 5), data = cur_data()),
      newdata = cur_data()
    )
  )

p111=ggplot(merged, aes(x = Year, y = fitted_GAM5_fs, color = subplot)) +
  geom_line(size = 1) +  # predicted trajectory
  geom_point(data = merged, 
             aes(x = Year, y = agb_gain, color = subplot),
             inherit.aes = FALSE, size = 2) +  # observed points
  theme_bw() +
  facet_wrap(~ subplot, scales = "free_y")+
  labs(x = "Year", 
       y = "agb_gain", title="GAM-k=5 fs") + theme(legend.position = "none")

p111 

plot(merged$agb_gain, merged$fitted_GAM5_fs)

p111.1 <- ggplot(merged, aes(Year, agb_gain - fitted_GAM5_fs)) +
  geom_point() +
  ggtitle("GAM residuals")

p111.1

ggplot(merged, aes(mean_z_tmax, agb_gain - fitted_GAM5_fs)) +
  geom_point() +
  geom_smooth(method = "lm", color = "red", se = FALSE) +
  facet_wrap(~ Treatment, scales = "free_y")+
  ggtitle("GAM5 bs = “fs” residuals vs Tmax")

#GAM captures early and late patterns faithfully, best edges behaviour. more variance on early years not linked to the model

slider

complete=TRUE (reject incomplete windows, NA at the beggining and the end), window per observation, not true “moving time windows”, window is just row-based.

merged <- merged %>%
  group_by(subplot) %>%
  arrange(Year) %>%
  mutate(fitted_slide_complete = slide_dbl(agb_gain, mean, .before = 2, .after = 2, .complete=TRUE))%>%
  ungroup()

p3=ggplot(merged, aes(x = Year, y = fitted_slide_complete, color = subplot)) +
  geom_line(size = 1) +  # predicted trajectory
  geom_point(data = merged, 
             aes(x = Year, y = agb_gain, color = subplot),
             inherit.aes = FALSE, size = 2) +  # observed points
  facet_wrap(~subplot, scales = "free_y")+
  theme_bw() +
  labs(x = "Year", 
       y = "agb_gain", title="slide - .before = 2, .after = 2, window per observation")+ theme(legend.position = "none")

p3 

plot(merged$agb_gain, merged$fitted_slide_complete)

p3.1 <- ggplot(merged, aes(Year, agb_gain - fitted_slide_complete)) +
  geom_point() +
  ggtitle("slide residuals")

p3.1

ggplot(merged, aes(mean_z_tmax, agb_gain - fitted_slide_complete)) +
  geom_point() +
  geom_smooth(method = "lm", color = "red", se = FALSE) +
  facet_wrap(~ Treatment, scales = "free_y")+
  ggtitle("slide_complete residuals vs Tmax")

complete=TRUE (reject incomplete windows, increased window NA at the beggining and the end), window per observation, not true “moving time windows”, window is just row-based.

merged <- merged %>%
  group_by(subplot) %>%
  arrange(Year) %>%
  mutate(fitted_slide_complete3 = slide_dbl(agb_gain, mean, .before = 3, .after = 3, .complete=TRUE))%>%
  ungroup()

p3=ggplot(merged, aes(x = Year, y = fitted_slide_complete3, color = subplot)) +
  geom_line(size = 1) +  # predicted trajectory
  geom_point(data = merged, 
             aes(x = Year, y = agb_gain, color = subplot),
             inherit.aes = FALSE, size = 2) +  # observed points
  facet_wrap(~subplot, scales = "free_y")+
  theme_bw() +
  labs(x = "Year", 
       y = "agb_gain", title="slide - .before = 3, .after = 3, window per observation")+ theme(legend.position = "none")

p3 

plot(merged$agb_gain, merged$fitted_slide_complete3)

p3.1 <- ggplot(merged, aes(Year, agb_gain - fitted_slide_complete3)) +
  geom_point() +
  ggtitle("slide residuals")

p3.1

ggplot(merged, aes(mean_z_tmax, agb_gain - fitted_slide_complete3)) +
  geom_point() +
  geom_smooth(method = "lm", color = "red", se = FALSE) +
  facet_wrap(~ Treatment, scales = "free_y")+
  ggtitle("slide_complete3 residuals vs Tmax")

**complete=TRUE, window per Year (NAs can be in the middle - anywhere Year gaps break the window)

merged <- merged %>%
  group_by(subplot) %>%
  arrange(Year) %>%
  mutate(fitted_slide_completeYear = slide_index_dbl(agb_gain, Year, mean, .before = 2, .after = 2, .complete=TRUE))%>%
  ungroup()

p33=ggplot(merged, aes(x = Year, y = fitted_slide_completeYear, color = subplot)) +
  geom_line(size = 1) +  # predicted trajectory
  geom_point(data = merged, 
             aes(x = Year, y = agb_gain, color = subplot),
             inherit.aes = FALSE, size = 2) +  # observed points
  facet_wrap(~subplot, scales= "free_y")+
  theme_bw() +
  labs(x = "Year", 
       y = "agb_gain", title="slide - .before = 2, .after = 2, window per year")+ theme(legend.position = "none")

p33 

plot(merged$agb_gain, merged$fitted_slide_completeYear)

p33.1 <- ggplot(merged, aes(Year, agb_gain - fitted_slide_completeYear)) +
  geom_point() +
  ggtitle("slide residuals")

p33.1

ggplot(merged, aes(mean_z_tmax, agb_gain - fitted_slide_completeYear)) +
  geom_point() +
  geom_smooth(method = "lm", color = "red", se = FALSE) +
  facet_wrap(~ Treatment, scales = "free_y")+
  ggtitle("slide_completeYear residuals vs Tmax")

**complete=TRUE, window per Year - increase window (more NAs)

merged <- merged %>%
  group_by(subplot) %>%
  arrange(Year) %>%
  mutate(fitted_slide_completeYear3 = slide_index_dbl(agb_gain, Year, mean, .before = 3, .after = 3, .complete=TRUE))%>%
  ungroup()

p33=ggplot(merged, aes(x = Year, y = fitted_slide_completeYear3, color = subplot)) +
  geom_line(size = 1) +  # predicted trajectory
  geom_point(data = merged, 
             aes(x = Year, y = agb_gain, color = subplot),
             inherit.aes = FALSE, size = 2) +  # observed points
  facet_wrap(~subplot, scales= "free_y")+
  theme_bw() +
  labs(x = "Year", 
       y = "agb_gain", title="slide - .before = 3, .after = 3, window per year")+ theme(legend.position = "none")

p33 

plot(merged$agb_gain, merged$fitted_slide_completeYear3)

p33.1 <- ggplot(merged, aes(Year, agb_gain - fitted_slide_completeYear3)) +
  geom_point() +
  ggtitle("slide residuals")

p33.1

ggplot(merged, aes(mean_z_tmax, agb_gain - fitted_slide_completeYear3)) +
  geom_point() +
  geom_smooth(method = "lm", color = "red", se = FALSE) +
  facet_wrap(~ Treatment, scales = "free_y")+
  ggtitle("slide_completeYear3 residuals vs Tmax")

complete=FALSE, no NAs, use less points

merged <- merged %>%
  group_by(subplot) %>%
  arrange(Year) %>%
  mutate(fitted_slide = slide_dbl(agb_gain, mean, .before = 2, .after = 2, .complete=FALSE))%>%
  ungroup()

p333=ggplot(merged, aes(x = Year, y = fitted_slide, color = subplot)) +
  geom_line(size = 1) +  # predicted trajectory
  geom_point(data = merged, 
             aes(x = Year, y = agb_gain, color = subplot),
             inherit.aes = FALSE, size = 2) +  # observed points
  facet_wrap(~subplot, scales = "free_y")+
  theme_bw() +
  labs(x = "Year", 
       y = "agb_gain", title="slide - .before = 2, .after = 2 - partial windows at edges")+ theme(legend.position = "none")

p333 

plot(merged$agb_gain, merged$fitted_slide)

p333.1 <- ggplot(merged, aes(Year, agb_gain - fitted_slide)) +
  geom_point() +
  ggtitle("slide residuals")


p333.1

ggplot(merged, aes(mean_z_tmax, agb_gain - fitted_slide)) +
  geom_point() +
  geom_smooth(method = "lm", color = "red", se = FALSE) +
  facet_wrap(~ Treatment, scales = "free_y")+
  ggtitle("slide residuals vs Tmax")

complete=FALSE, increased window no NAs, use less points

merged <- merged %>%
  group_by(subplot) %>%
  arrange(Year) %>%
  mutate(fitted_slide3 = slide_dbl(agb_growth, mean, .before = 3, .after = 3, .complete=FALSE))%>%
  ungroup()

p333=ggplot(merged, aes(x = Year, y = fitted_slide3, color = subplot)) +
  geom_line(size = 1) +  # predicted trajectory
  geom_point(data = merged, 
             aes(x = Year, y = agb_growth, color = subplot),
             inherit.aes = FALSE, size = 2) +  # observed points
  facet_wrap(~subplot, scales = "free_y")+
  theme_bw() +
  labs(x = "Year", 
       y = "agb_growth", title="slide - .before = 3, .after = 3 - partial windows at edges")+ theme(legend.position = "none")

p333 

plot(merged$agb_gain, merged$fitted_slide3)

p333.1 <- ggplot(merged, aes(Year, agb_gain - fitted_slide3)) +
  geom_point() +
  ggtitle("slide residuals")


p333.1

ggplot(merged, aes(mean_z_tmax, agb_gain - fitted_slide3)) +
  geom_point() +
  geom_smooth(method = "lm", color = "red", se = FALSE) +
  facet_wrap(~ Treatment, scales = "free_y")+
  ggtitle("slide3 residuals vs Tmax")

root mean square error

rmse_results <- merged %>%
  group_by(subplot) %>%
  summarise(
    RMSE_GAM8   = sqrt(mean((agb_gain - fitted_GAM8)^2, na.rm = TRUE)),
    RMSE_GAM5_fs   = sqrt(mean((agb_gain - fitted_GAM5_fs)^2, na.rm = TRUE)),
    RMSE_GAM5   = sqrt(mean((agb_gain - fitted_GAM5)^2, na.rm = TRUE)),
    RMSE_SLIDE  = sqrt(mean((agb_gain - fitted_slide)^2, na.rm = TRUE)),
    RMSE_SLIDE_COMPLETE  = sqrt(mean((agb_gain - fitted_slide_complete)^2, na.rm = TRUE)),
    RMSE_SLIDE_COMPLETE_YEAR  = sqrt(mean((agb_gain - fitted_slide_completeYear)^2, na.rm = TRUE)),
    RMSE_SLIDE3  = sqrt(mean((agb_gain - fitted_slide3)^2, na.rm = TRUE)),
    RMSE_SLIDE_COMPLETE3  = sqrt(mean((agb_gain - fitted_slide_complete3)^2, na.rm = TRUE)),
    RMSE_SLIDE_COMPLETE_YEAR3  = sqrt(mean((agb_gain - fitted_slide_completeYear3)^2, na.rm = TRUE)),
  ) %>% 
  ungroup()

rmse_results <- rmse_results %>%
  mutate(
    best_method = colnames(.)[
      apply(select(., -subplot), 1, function(x) which.min(replace(x, is.na(x), Inf))) + 1
    ]
  )

as.data.table(rmse_results)

    subplot RMSE_GAM8 RMSE_GAM5_fs RMSE_GAM5 RMSE_SLIDE RMSE_SLIDE_COMPLETE
     <fctr>     <num>        <num>     <num>      <num>               <num>
 1:     1_1 0.8242306    0.8533249 0.8242306  0.8512805           0.7378002
 2:     1_2 0.5131833    0.5789256 0.5131833  0.5372679           0.5699395
 3:     1_3 0.7635309    0.7874809 0.7635309  0.8191958           0.7836013
 4:     1_4 0.6408596    0.6424491 0.6408596  0.6790087           0.6760744
 5:    10_1 0.9446733    1.0575580 0.8989452  0.9372639           0.9334989
 6:    10_2 0.8962299    0.9570244 0.8949176  0.9131286           0.8685274
 7:    10_3 0.8212653    0.9543175 0.8591196  0.8285771           0.7875102
 8:    10_4 0.7843370    0.9680006 0.7968818  0.8269127           0.8491458
 9:    11_1 0.7517646    0.7989772 0.7512557  0.7894937           0.6773169
10:    11_2 0.7420913    0.7455335 0.7420913  0.7192910           0.6923522
11:    11_3 0.7985344    0.7744485 0.7985344  0.8088344           0.6277109
12:    11_4 0.7648175    0.7848663 0.7647166  0.6946322           0.5357400
13:    12_1 1.0000814    1.0705614 0.9923200  0.8876736           0.8323566
14:    12_2 0.8279142    0.8478490 0.8279142  0.8476785           0.7999676
15:    12_3 1.0552169    1.1000136 1.0411502  0.9760758           0.9699883
16:    12_4 0.9879689    0.9707440 0.9879689  0.9890025           0.9754880
17:    13_1 1.1053452    1.1165349 1.1111648  0.8691264           0.8419337
18:    13_2 0.9532163    0.9420210 0.9532163  0.9443286           0.9011302
19:    13_3 0.6815448    0.8210827 0.6878990  0.7330095           0.4256024
20:    13_4 0.8537066    0.9028229 0.8477012  0.7674951           0.5731913
21:    14_1 0.8211036    0.8370600 0.8227841  0.8258056           0.7549048
22:    14_2 0.4955378    0.5181103 0.4976393  0.4754689           0.4167890
23:    14_3 0.6775247    0.7076194 0.6775874  0.6120385           0.5713863
24:    14_4 0.6810291    0.6683124 0.6809399  0.5504938           0.4500254
25:    15_1 0.8015749    0.8138566 0.8015749  0.6479123           0.5671449
26:    15_2 0.7226396    0.7130944 0.7226396  0.6633465           0.7088132
27:    15_3 0.7490577    0.7351057 0.7490577  0.7687220           0.6902662
28:    15_4 1.0153011    1.0026657 1.0163212  0.8211689           0.7807313
29:     2_1 0.3719264    0.6628040 0.6543560  0.5618052           0.4880375
30:     2_2 0.3677282    0.5109481 0.5237894  0.4514330           0.4861201
31:     2_3 0.7239030    0.8183051 0.7239030  0.7500960           0.6806372
32:     2_4 0.6862915    0.6745867 0.6862915  0.6405031           0.6964029
33:     3_1 0.7758756    0.7622791 0.7758756  0.6897714           0.7481650
34:     3_2 0.4334168    0.5214514 0.5026218  0.4658375           0.4093548
35:     3_3 0.7653404    0.8417284 0.7653404  0.7108028           0.7554039
36:     3_4 0.7379422    0.7191307 0.7379422  0.6348814           0.6412186
37:     4_1 0.3655344    0.6028095 0.3660401  0.3726241           0.3957978
38:     4_2 0.5230318    0.5797229 0.5092958  0.5505045           0.6106932
39:     4_3 0.7726508    1.0006275 0.8001743  0.7703836           0.6512226
40:     4_4 0.5493570    0.7349730 0.5554626  0.5859284           0.6249544
41:     5_1 0.5840682    0.7008998 0.6721414  0.5678325           0.6063996
42:     5_2 0.5595106    0.6214271 0.5636336  0.5066311           0.5672482
43:     5_3 0.8354678    0.8292434 0.8354678  0.7972868           0.7720176
44:     5_4 0.8459623    0.8553703 0.8459623  0.8029789           0.9043768
45:     6_1 1.0292794    1.0182638 1.0306605  1.0025727           0.9289293
46:     6_2 0.6629219    0.8794186 0.8837724  0.8158759           0.7809733
47:     6_3 0.9188234    0.9395424 0.9188234  0.9411167           0.6718756
48:     6_4 0.6166674    0.8674358 0.8281162  0.7823398           0.6733959
49:     7_1 0.7178275    0.7089129 0.7178275  0.6548787           0.6880387
50:     7_2 0.7422772    0.7377438 0.7422772  0.7933843           0.7995247
51:     7_3 0.6555534    0.6637050 0.6555534  0.6681872           0.6957817
52:     7_4 0.6897942    0.6865442 0.6897942  0.6382536           0.6589007
53:     8_1 0.5606565    0.7403506 0.5528924  0.5051216           0.5072790
54:     8_2 0.9719646    0.9744401 0.9719646  0.9026908           0.9937535
55:     8_3 0.6554025    0.7395910 0.6499500  0.6316303           0.6575989
56:     8_4 0.7739914    0.9272844 0.7610548  0.8138452           0.8451426
57:     9_1 0.5843730    0.5890522 0.5843730  0.5627049           0.5703913
58:     9_2 0.8686738    0.8580513 0.8682625  0.8583545           0.8803016
59:     9_3 0.6991372    0.6842350 0.6993837  0.6661502           0.6063643
60:     9_4 0.7440635    0.7378926 0.7440635  0.7597730           0.7076661
    subplot RMSE_GAM8 RMSE_GAM5_fs RMSE_GAM5 RMSE_SLIDE RMSE_SLIDE_COMPLETE
    RMSE_SLIDE_COMPLETE_YEAR RMSE_SLIDE3 RMSE_SLIDE_COMPLETE3
                       <num>       <num>                <num>
 1:                0.7903062   0.8631095            0.7052939
 2:                0.4786826   0.6275165            0.5758997
 3:                0.6924043   0.8348774            0.6791101
 4:                0.6510481   0.7076796            0.7117678
 5:                0.9958413   1.5384507            0.9609460
 6:                0.9109461   1.2344246            0.7445580
 7:                0.8674531   1.2156838            0.9184014
 8:                0.8232213   1.1934298            0.8495921
 9:                0.7459328   0.7753979            0.6509495
10:                0.7363253   0.7673774            0.6679234
11:                0.6705234   0.7789627            0.6446485
12:                0.5999013   0.8245071            0.5828276
13:                1.0790364   1.4179509            1.0766924
14:                0.8622936   1.1232582            0.8799680
15:                1.1530924   1.5519536            1.0277964
16:                1.0468404   1.3183472            1.0463501
17:                0.9784808   1.0352346            1.1369303
18:                0.9808817   0.9912271            0.9793851
19:                0.6625227   0.8938999            0.5540179
20:                0.7670313   0.9538784            0.7039857
21:                0.8135286   0.8767844            0.7675604
22:                0.3133097   0.5329811            0.3996454
23:                0.5981214   0.6083958            0.5064219
24:                0.5543879   0.6378018            0.4597678
25:                0.5857080   0.6645573            0.5640312
26:                0.5662909   0.7009633            0.6565938
27:                0.6689654   0.8206176            0.7306196
28:                0.8403002   0.9788645            0.9382542
29:                0.3698907   0.9578987            0.5985820
30:                0.3780277   0.7456711            0.5505659
31:                0.6541470   1.1167539            0.5351314
32:                0.7294166   0.8302098            0.6680572
33:                0.7695540   1.3533010            0.6926796
34:                0.3210781   1.0612416            0.4090400
35:                0.7578459   1.5869044            0.6464870
36:                0.7986745   1.3256850            0.6537454
37:                0.4308294   1.0110476            0.4259681
38:                0.5663120   1.2699084            0.5480453
39:                0.8858252   1.6454824            0.7343123
40:                0.5841109   1.2132814            0.6074167
41:                0.5804404   1.1204944            0.7062864
42:                0.5464812   0.8464098            0.6241970
43:                0.8536224   1.0267702            0.8187822
44:                0.9302481   1.0923013            0.9201969
45:                0.8407195   1.0914316            0.9528914
46:                0.6695340   0.9415860            0.7881393
47:                0.8098509   0.9499470            0.6771344
48:                0.6369743   0.9683585            0.7234772
49:                0.6666748   0.9084604            0.7334917
50:                0.7283272   0.9100858            0.6968475
51:                0.7309603   0.8203331            0.7408434
52:                0.8295693   0.8178553            0.8037447
53:                0.6971338   1.4236395            0.7243757
54:                1.0906926   1.9622760            1.1508375
55:                0.7453144   1.3036657            0.5974037
56:                0.7455778   1.6741037            0.8090783
57:                0.5562750   0.9102644            0.5981279
58:                0.9136981   1.0961633            0.9055162
59:                0.6285110   0.8772287            0.6498184
60:                0.6581235   1.0790123            0.7228062
    RMSE_SLIDE_COMPLETE_YEAR RMSE_SLIDE3 RMSE_SLIDE_COMPLETE3
    RMSE_SLIDE_COMPLETE_YEAR3               best_method
                        <num>                    <char>
 1:                 0.7363183      RMSE_SLIDE_COMPLETE3
 2:                 0.4838016  RMSE_SLIDE_COMPLETE_YEAR
 3:                 0.7241385      RMSE_SLIDE_COMPLETE3
 4:                 0.6349536 RMSE_SLIDE_COMPLETE_YEAR3
 5:                 1.0502722                 RMSE_GAM5
 6:                 0.8217983      RMSE_SLIDE_COMPLETE3
 7:                 0.9480886       RMSE_SLIDE_COMPLETE
 8:                 0.8929905                 RMSE_GAM8
 9:                 0.7178193      RMSE_SLIDE_COMPLETE3
10:                 0.7620733      RMSE_SLIDE_COMPLETE3
11:                 0.6669592       RMSE_SLIDE_COMPLETE
12:                 0.6537341       RMSE_SLIDE_COMPLETE
13:                 1.1209510       RMSE_SLIDE_COMPLETE
14:                 0.9039531       RMSE_SLIDE_COMPLETE
15:                 1.2161317       RMSE_SLIDE_COMPLETE
16:                 1.1549717              RMSE_GAM5_fs
17:                 1.0858397       RMSE_SLIDE_COMPLETE
18:                 0.9670784       RMSE_SLIDE_COMPLETE
19:                 0.5125977       RMSE_SLIDE_COMPLETE
20:                 0.7386343       RMSE_SLIDE_COMPLETE
21:                 0.8607089       RMSE_SLIDE_COMPLETE
22:                 0.3541762  RMSE_SLIDE_COMPLETE_YEAR
23:                 0.5230778      RMSE_SLIDE_COMPLETE3
24:                 0.4972726       RMSE_SLIDE_COMPLETE
25:                 0.5802532      RMSE_SLIDE_COMPLETE3
26:                 0.4630949 RMSE_SLIDE_COMPLETE_YEAR3
27:                 0.7089224  RMSE_SLIDE_COMPLETE_YEAR
28:                 0.9225724       RMSE_SLIDE_COMPLETE
29:                 0.3311902 RMSE_SLIDE_COMPLETE_YEAR3
30:                 0.4431806                 RMSE_GAM8
31:                 0.4325432 RMSE_SLIDE_COMPLETE_YEAR3
32:                 0.7118625                RMSE_SLIDE
33:                 0.7752400                RMSE_SLIDE
34:                 0.3183548 RMSE_SLIDE_COMPLETE_YEAR3
35:                 0.6782628      RMSE_SLIDE_COMPLETE3
36:                 0.7211694                RMSE_SLIDE
37:                 0.4227364                 RMSE_GAM8
38:                 0.5634202                 RMSE_GAM5
39:                 0.7684061       RMSE_SLIDE_COMPLETE
40:                 0.6262974                 RMSE_GAM8
41:                 0.6219605                RMSE_SLIDE
42:                 0.5611753                RMSE_SLIDE
43:                 0.7614985 RMSE_SLIDE_COMPLETE_YEAR3
44:                 0.9638445                RMSE_SLIDE
45:                 0.9146083  RMSE_SLIDE_COMPLETE_YEAR
46:                 0.7639693                 RMSE_GAM8
47:                 0.6325150 RMSE_SLIDE_COMPLETE_YEAR3
48:                 0.7192003                 RMSE_GAM8
49:                 0.6324573 RMSE_SLIDE_COMPLETE_YEAR3
50:                 0.6610246 RMSE_SLIDE_COMPLETE_YEAR3
51:                 0.7948062                 RMSE_GAM8
52:                 0.8676862                RMSE_SLIDE
53:                 0.7443500                RMSE_SLIDE
54:                 1.2534069                RMSE_SLIDE
55:                 0.7480171      RMSE_SLIDE_COMPLETE3
56:                 0.8031148  RMSE_SLIDE_COMPLETE_YEAR
57:                 0.5309625 RMSE_SLIDE_COMPLETE_YEAR3
58:                 0.8576640 RMSE_SLIDE_COMPLETE_YEAR3
59:                 0.5542777 RMSE_SLIDE_COMPLETE_YEAR3
60:                 0.6452868 RMSE_SLIDE_COMPLETE_YEAR3
    RMSE_SLIDE_COMPLETE_YEAR3               best_method

lm

vars <- grep("^fitted_", names(merged), value = TRUE)

results <- data.frame(
  method = vars,
  intercept = NA,
  slope = NA,
  p_value = NA
)



for (i in seq_along(vars)) {
  
  v <- vars[i]
  model <- lm(as.formula(paste0(v, " - agb_gain ~ Year")), data = merged)
  s <- summary(model)
  
  results$intercept[i] <- coef(s)[1, 1]
  results$slope[i]     <- coef(s)[2, 1]
  results$p_value[i]   <- coef(s)[2, 4]
}

results

                      method     intercept         slope      p_value
1                fitted_GAM5  4.303302e-11 -2.146657e-14 1.000000e+00
2                fitted_GAM8 -5.712758e-12  5.704802e-15 1.000000e+00
3             fitted_GAM5_fs -2.024540e-01  1.009893e-04 9.690271e-01
4      fitted_slide_complete  9.869937e+00 -4.919834e-03 1.265646e-01
5     fitted_slide_complete3  1.547904e+01 -7.730575e-03 6.468196e-02
6  fitted_slide_completeYear -5.276129e+00  2.655359e-03 3.813579e-01
7 fitted_slide_completeYear3  1.407310e+01 -7.015788e-03 6.100874e-02
8               fitted_slide  2.711661e+00 -1.354328e-03 5.766612e-01
9              fitted_slide3 -4.654542e+01  2.290721e-02 5.291172e-15

temporal signal on residuals

methods <- c(
  "fitted_GAM5",
  "fitted_GAM8",
  "fitted_GAM5_fs",
  "fitted_slide_completeYear",
  "fitted_slide_complete",
  "fitted_slide",
  "fitted_slide3",
  "fitted_slide_completeYear3",
  "fitted_slide_complete3"
)

for (m in methods) {
  merged[[paste0("resid_", m)]] <- merged$agb_gain- merged[[m]]
}

pairwise_corr <- function(df, col) {
  wide <- tidyr::pivot_wider(df[, c("subplot","Year", col)],
                             names_from = subplot, values_from = all_of(col))
  M <- as.matrix(wide[, -1])          # remove Year
  C <- cor(M, use = "pairwise.complete.obs")
  mean(C[upper.tri(C)], na.rm = TRUE)
}

temporal_signal <- data.frame(
  method = methods,
  mean_pairwise_corr = sapply(paste0("resid_", methods),
                              function(x) pairwise_corr(merged, x))
)

temporal_signal

                                                     method mean_pairwise_corr
resid_fitted_GAM5                               fitted_GAM5          0.3236167
resid_fitted_GAM8                               fitted_GAM8          0.3256980
resid_fitted_GAM5_fs                         fitted_GAM5_fs          0.3073723
resid_fitted_slide_completeYear   fitted_slide_completeYear          0.2552291
resid_fitted_slide_complete           fitted_slide_complete          0.2806969
resid_fitted_slide                             fitted_slide          0.3243473
resid_fitted_slide3                           fitted_slide3          0.3258910
resid_fitted_slide_completeYear3 fitted_slide_completeYear3          0.2569271
resid_fitted_slide_complete3         fitted_slide_complete3          0.3262622

save(merged, file="allplots_smooth_models_agb_gain.rda")